function [pigrid,ssrlog,survivorfrac] = ASPS(pigrid0,knotvec_c,P,BSstate,keepnum,partial)
% % % % % ASPS STANDS FOR "ADAPTIVE START-POINT SEARCH"


% test=[]; scalefactor = 0.5; stepsize = 4/5 ; steps = 12 ; for ii=1:steps; test = [test scalefactor*(stepsize^(ii-1))]; end;
% scalefactor = 0.375; stepsize = 3/4 ; steps = 10 ; for ii=1:steps; test = [test scalefactor*(stepsize^(ii-1))]; end;
% scalefactor = 0.25; stepsize = 2/3 ; steps = 8 ; for ii=1:steps; test = [test scalefactor*(stepsize^(ii-1))]; end;
% scalefactor = 0.125; stepsize = 1/2 ; steps = 6 ; for ii=1:steps; test = [test scalefactor*(stepsize^(ii-1))]; end;
% scalefactor = 0.0625; stepsize = 1/4 ; steps = 4 ; for ii=1:steps; test = [test scalefactor*(stepsize^(ii-1))]; end;
% test=sort(1-test); scatter(test,ones(size(test)),'filled'); grid on; box on; disp(test'); disp(size(test))


if nargin<6
    partial = 0 ;
end
if BSstate(1)==0
    feastol = 10^-8 ;
else
    feastol = 10^-8 ;
end


student_id   = P.student_id ; 
Q            = P.Q ; 
T            = P.T ; 
ThetaE       = P.ThetaE ; 
contract     = P.contract ; 
Q80freq      = P.Q80freq ; 
pay          = P.pay ; 
SimTimeCum   = P.SimTimeCum ; 
Activeids    = P.Activeids ; 
Marginalids1 = P.Marginalids1 ; 
Marginalids2 = P.Marginalids2 ; 
Tdom         = P.Tdom ; 
EQgivenT_lb  = P.EQgivenT_lb ; 
EQgivenT_ub  = P.EQgivenT_ub ; 
cutoffs      = P.cutoffs ; 
Aeq          = P.Aeq ; 
beq          = P.beq ; 
Aineq        = P.Aineq ; 
bineq        = P.bineq ; 
lb           = P.lb ; 
ub           = P.ub ; 

%%%%The next few lines check to make sure that the pigrid0 matrix has the correct format:
[testrow,testcol] = size(pigrid0) ;
if testcol==length(knotvec_c)+2
    pigrid0 = pigrid0' ;
elseif testrow~=length(knotvec_c)+2
    error('ASPS:BadInput','pigrid0 SHOULD BE COMPRISED OF PARAMETER VECTORS OF LENGTH length(knotvec_c)+2...')
end



if BSstate(1)==0 || BSstate(1)==10 || BSstate(1)==11  
    betaq1    = P.betaq1 ; 
    betaq2    = P.betaq2 ; 
    betaq3    = P.betaq3 ; 
    knotvec_q = P.knotvec_q ; 
else
    betaq1    = [] ; 
    betaq2    = [] ; 
    betaq3    = [] ; 
    knotvec_q = [] ; 
end

disp('%%%%000000000000000000000000000000000000%%%%')
disp('%%%% ADAPTIVE START POINT SEARCH, ROUND 0...')
disp('%%%%000000000000000000000000000000000000%%%%')
    %%%%Weed out any repeats:
    pigrid0  = unique(pigrid0','rows','stable')' ;
    %%%%Weed out any vectors that violat feasibility:
    [eqindx,ineqindx,boundindx]  = feasiblecheck(pigrid0,Aeq,beq,Aineq,bineq,lb,ub,feastol) ; 
    pigrid0(1,eqindx==0) = beq(1) ; 
    pigrid0(2,eqindx==0) = beq(2) ; 
    if sum(ineqindx==0 | boundindx==0)>0
        pigrid0 = pigrid0(:,ineqindx==1 & boundindx==1) ;
        disp('*********************************************************************************************')
        disp('*********************************************************************************************')
        disp(strcat('* THERE WERE_',num2str(sum(ineqindx==0)),'_INFEASIBLE INPUT VECTORS DISCARDED, BASED ON INEQUALITY/BOUND CONSTRAINTS... *'))  
        disp('*********************************************************************************************')
        disp('*********************************************************************************************')
    end 
    pi_c     = pigrid0(:,1) ;  %%%%THIS LINE IS JUST FOR DISPLAY PURPOSES, BUT THE FIRST COLUMN OF pigrid SHOULD CONTAIN THE CURRENT BEST GUESS...
    ssrlog0  = nan(length(pigrid0(1,:)),1) ;
    disp(strcat('evaluating_',num2str(length(ssrlog0)),'_candidate parameter vectors.'))
    tic  
    parfor ii=1:length(pigrid0(1,:))
        pitemp  = pigrid0(:,ii) ;
        ssrtemp = ObjFun_cMSM(pitemp,knotvec_c,student_id,Q,T,...  %%%
                               ThetaE,contract,Q80freq,pay,SimTimeCum,... %%%
                               Activeids,Marginalids1,Marginalids2,...
                               Tdom,EQgivenT_lb,EQgivenT_ub,cutoffs,... %%%
                               betaq1,betaq2,betaq3,knotvec_q,0,BSstate) ;
        ssrlog0(ii)  = ssrtemp ;
    end
    toc  
    nextindx = find(ssrlog0==min(ssrlog0),1,'first') ; 
    if length(nextindx)>1
        if max(abs(pigrid0(:,nextindx(1))-pigrid0(:,nextindx(2)))./pigrid0(:,nextindx(1)))<10^-7
            pinext = pigrid0(:,nextindx(1)) ;
        else
            warning('Two minimizers found where relative difference is less than 10^-7...')
            pinext = pigrid0(:,nextindx(1)) ;
        end
    else
        pinext   = pigrid0(:,nextindx) ;
    end
    format shortg
    CheckAltGuesses = [pi_c pinext 100*(pinext-pi_c)./pi_c; ssrlog0(1) ssrlog0(nextindx(1)) 100*(ssrlog0(nextindx(1))-ssrlog0(1))/ssrlog0(1)]   %#ok
    pi_c   = pinext ; 
    pigrid = pigrid0 ; 
    ssrlog = ssrlog0 ; 

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%  THE FOLLOWING BLOCK OF CODE INTEGRATES AN ADAPTIVE GRID SEARCH ROUTINE WITH THE GMM
    %%%%%%%%%%  SOLVER TO FIND A SUITABLE INITIAL GUESS AND THEN DOUBLE CHECK TO SEE IF THE SOLUTION
    %%%%%%%%%%  THE SOLVER CONVERGES TO CAN BE IMPROVED UPON.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    startvaldisplay = 0 ;
    disp('%%%%111111111111111111111111111111111111%%%%')
    disp('%%%% ADAPTIVE START POINT SEARCH, ROUND 1...')
    disp('%%%%111111111111111111111111111111111111%%%%')
    scalefactor     = 0.5 ;
    upperlinwt      = 1 ;
    steps           =5 ;
    stepsize        = 4/5 ;
    if length(knotvec_c)>=10
        uppercurvfactor = 0.2 ;
        pigrid_supp = startvalgeneratorMultiStep(knotvec_c,pi_c,scalefactor,upperlinwt,uppercurvfactor,steps,stepsize,4,startvaldisplay) ;
    elseif length(knotvec_c)==9
        uppercurvfactor = 0.3 ;
        pigrid_supp = startvalgeneratorMultiStep(knotvec_c,pi_c,scalefactor,upperlinwt,uppercurvfactor,steps,stepsize,3,startvaldisplay) ;
    elseif length(knotvec_c)==8
        uppercurvfactor = 0.4 ;
        pigrid_supp = startvalgeneratorMultiStep(knotvec_c,pi_c,scalefactor,upperlinwt,uppercurvfactor,steps,stepsize,2,startvaldisplay) ;
    end
    %%%%Weed out any repeats:
    pigrid_supp = unique(pigrid_supp','rows','stable')' ;
    %%%%Weed out any vectors that violate feasibility:
    [eqindx,ineqindx,boundindx]  = feasiblecheck(pigrid_supp,Aeq,beq,Aineq,bineq,lb,ub,feastol) ; 
    pigrid_supp(1,eqindx==0) = beq(1) ; 
    pigrid_supp(2,eqindx==0) = beq(2) ; 
    if sum(ineqindx==0|boundindx==0)>0 
        pigrid_supp = pigrid_supp(:,ineqindx==1&boundindx==1) ;
    end 
    ssrlog_supp  = nan(length(pigrid_supp(1,:)),1) ;
    disp(strcat('evaluating_',num2str(length(ssrlog_supp)),'_candidate parameter vectors.')) 
    tic
    parfor ii=1:length(pigrid_supp(1,:)) 
        pitemp = pigrid_supp(:,ii) ;
        ssrtemp = ObjFun_cMSM(pitemp,knotvec_c,student_id,Q,T,...  %%%
                               ThetaE,contract,Q80freq,pay,SimTimeCum,... %%%
                               Activeids,Marginalids1,Marginalids2,...
                               Tdom,EQgivenT_lb,EQgivenT_ub,cutoffs,... %%%
                               betaq1,betaq2,betaq3,knotvec_q,0,BSstate) ;
        ssrlog_supp(ii)  = ssrtemp ;
    end
    toc
    pigrid(:,length(pigrid(1,:))+1:length(pigrid(1,:))+length(pigrid_supp(1,:))) = pigrid_supp ;  
    ssrlog(length(ssrlog)+1:length(ssrlog)+length(ssrlog_supp))                  = ssrlog_supp ;
    nextindx = find(ssrlog==min(ssrlog),1,'first') ;
    if length(nextindx)>1
        if max(abs(pigrid(:,nextindx(1))-pigrid(:,nextindx(2)))./pigrid(:,nextindx(1)))<10^-7
            pinext = pigrid(:,nextindx(1)) ;
        else
            warning('Two minimizers found where relative difference is less than 10^-7...')
            pinext = pigrid(:,nextindx(1)) ;
        end
    else
        pinext   = pigrid(:,nextindx) ;
    end
    format shortg
    CheckAltGuesses = [pi_c pinext 100*(pinext-pi_c)./pi_c; ssrlog(1) ssrlog(nextindx(1)) 100*(ssrlog(nextindx(1))-ssrlog(1))/ssrlog(1)]   %#ok
    pi_c = pinext ;
    if partial==0
        %%%%Now, we re-do the grid search using the best guess from the previous round and more steps
        disp('%%%%222222222222222222222222222222222222%%%%')
        disp('%%%% ADAPTIVE START POINT SEARCH, ROUND 2...')
        disp('%%%%222222222222222222222222222222222222%%%%')
        scalefactor     = 0.375 ;
        upperlinwt      = 0.75 ;
        uppercurvfactor = uppercurvfactor*(7/8) ;
        steps           = 5 ;
        stepsize        = 3/4 ;
        if length(knotvec_c)>=10
            pigrid_supp   = startvalgeneratorMultiStep(knotvec_c,pi_c,scalefactor,upperlinwt,uppercurvfactor,steps,stepsize,4,startvaldisplay) ;
        elseif length(knotvec_c)==9
            pigrid_supp   = startvalgeneratorMultiStep(knotvec_c,pi_c,scalefactor,upperlinwt,uppercurvfactor,steps,stepsize,3,startvaldisplay) ;
        elseif length(knotvec_c)==8
            pigrid_supp   = startvalgeneratorMultiStep(knotvec_c,pi_c,scalefactor,upperlinwt,uppercurvfactor,steps,stepsize,2,startvaldisplay) ;
        end
        %%%%Weed out any repeats:
        pigrid_supp = unique(pigrid_supp','rows','stable')' ;
        %%%%Weed out any vectors that violate feasibility:
        [eqindx,ineqindx,boundindx]  = feasiblecheck(pigrid_supp,Aeq,beq,Aineq,bineq,lb,ub,feastol) ;
        pigrid_supp(1,eqindx==0) = beq(1) ;
        pigrid_supp(2,eqindx==0) = beq(2) ;
        if sum(ineqindx==0|boundindx==0)>0
            pigrid_supp = pigrid_supp(:,ineqindx==1&boundindx==1) ;
        end
        ssrlog_supp   = nan(length(pigrid_supp(1,:)),1) ;
        disp(strcat('evaluating_',num2str(length(ssrlog_supp)),'_candidate parameter vectors.'))
        tic
        parfor ii=1:length(pigrid_supp(1,:))
            try
                pitemp = pigrid_supp(:,ii) ;
                ssrtemp = ObjFun_cMSM(pitemp,knotvec_c,student_id,Q,T,...  %%%
                    ThetaE,contract,Q80freq,pay,SimTimeCum,... %%%
                    Activeids,Marginalids1,Marginalids2,...
                    Tdom,EQgivenT_lb,EQgivenT_ub,cutoffs,... %%%
                    betaq1,betaq2,betaq3,knotvec_q,0,BSstate) ;
                ssrlog_supp(ii)  = ssrtemp ;
            catch
                ssrlog_supp(ii) = inf ;
            end
        end
        toc
        pigrid(:,length(pigrid(1,:))+1:length(pigrid(1,:))+length(pigrid_supp(1,:))) = pigrid_supp ;
        ssrlog(length(ssrlog)+1:length(ssrlog)+length(ssrlog_supp))                  = ssrlog_supp ;
        nextindx = find(ssrlog==min(ssrlog),1,'first') ;
        if length(nextindx)>1
            if max(abs(pigrid(:,nextindx(1))-pigrid(:,nextindx(2)))./pigrid(:,nextindx(1)))<10^-7
                pinext = pigrid(:,nextindx(1)) ;
            else
                warning('Two minimizers found where relative difference is less than 10^-7...')
                pinext = pigrid(:,nextindx(1)) ;
            end
        else
            pinext   = pigrid(:,nextindx) ;
        end
        format shortg
        CheckAltGuesses = [pi_c pinext 100*(pinext-pi_c)./pi_c; ssrlog(1) ssrlog(nextindx(1)) 100*(ssrlog(nextindx(1))-ssrlog(1))/ssrlog(1)]   %#ok
        pi_c = pinext ;
    end
    %%%%Now, we re-do the grid search once using the best guess from the previous round, a smaller
    %%%%stepsize, and smaller start values:
    disp('%%%%333333333333333333333333333333333333%%%%')
    disp('%%%% ADAPTIVE START POINT SEARCH, ROUND 3...')
    disp('%%%%333333333333333333333333333333333333%%%%')
    scalefactor     = 0.25 ;
    upperlinwt      = 0.5 ;
    uppercurvfactor = uppercurvfactor*(6/7) ; %%%%This implies (3/4) of the original value
    steps           = 3 ;
    stepsize        = 2/3 ;
    if length(knotvec_c)>=10
        pigrid_supp   = startvalgeneratorMultiStep(knotvec_c,pi_c,scalefactor,upperlinwt,uppercurvfactor,steps,stepsize,4,startvaldisplay) ;
    elseif length(knotvec_c)==9
        pigrid_supp   = startvalgeneratorMultiStep(knotvec_c,pi_c,scalefactor,upperlinwt,uppercurvfactor,steps,stepsize,3,startvaldisplay) ;
    elseif length(knotvec_c)==8
        pigrid_supp   = startvalgeneratorMultiStep(knotvec_c,pi_c,scalefactor,upperlinwt,uppercurvfactor,steps,stepsize,2,startvaldisplay) ;
    end
    %%%%Weed out any repeats:
    pigrid_supp = unique(pigrid_supp','rows','stable')' ;
    %%%%Weed out any vectors that violate feasibility:
    [eqindx,ineqindx,boundindx]  = feasiblecheck(pigrid_supp,Aeq,beq,Aineq,bineq,lb,ub,feastol) ; 
    pigrid_supp(1,eqindx==0) = beq(1) ; 
    pigrid_supp(2,eqindx==0) = beq(2) ; 
    if sum(ineqindx==0|boundindx==0)>0 
        pigrid_supp = pigrid_supp(:,ineqindx==1&boundindx==1) ;
    end
    ssrlog_supp   = nan(length(pigrid_supp(1,:)),1) ;
    disp(strcat('evaluating_',num2str(length(ssrlog_supp)),'_candidate parameter vectors.')) 
    tic
    parfor ii=1:length(pigrid_supp(1,:))
        try
            pitemp = pigrid_supp(:,ii) ;
            ssrtemp = ObjFun_cMSM(pitemp,knotvec_c,student_id,Q,T,...  %%%
                                   ThetaE,contract,Q80freq,pay,SimTimeCum,... %%%
                                   Activeids,Marginalids1,Marginalids2,...
                                   Tdom,EQgivenT_lb,EQgivenT_ub,cutoffs,... %%%
                                   betaq1,betaq2,betaq3,knotvec_q,0,BSstate) ;
            ssrlog_supp(ii)  = ssrtemp ;
        catch
            ssrlog_supp(ii) = inf ;
        end
    end
    toc
    pigrid(:,length(pigrid(1,:))+1:length(pigrid(1,:))+length(pigrid_supp(1,:))) = pigrid_supp ;  
    ssrlog(length(ssrlog)+1:length(ssrlog)+length(ssrlog_supp))                  = ssrlog_supp ;
    nextindx = find(ssrlog==min(ssrlog),1,'first') ;
    if length(nextindx)>1
        if max(abs(pigrid(:,nextindx(1))-pigrid(:,nextindx(2)))./pigrid(:,nextindx(1)))<10^-7
            pinext = pigrid(:,nextindx(1)) ;
        else
            warning('Two minimizers found where relative difference is less than 10^-7...')
            pinext = pigrid(:,nextindx(1)) ;
        end
    else
        pinext   = pigrid(:,nextindx) ;
    end
    format shortg
    CheckAltGuesses = [pi_c pinext 100*(pinext-pi_c)./pi_c; ssrlog(1) ssrlog(nextindx(1)) 100*(ssrlog(nextindx(1))-ssrlog(1))/ssrlog(1)]   %#ok
    pi_c = pinext ;
    if partial==0
        %%%%Now, we re-do the grid search once using the best guess from the previous round, a smaller
        %%%%stepsize, and smaller start values:
        disp('%%%%444444444444444444444444444444444444%%%%')
        disp('%%%% ADAPTIVE START POINT SEARCH, ROUND 4...')
        disp('%%%%444444444444444444444444444444444444%%%%')
        scalefactor     = 0.125 ;
        upperlinwt      = upperlinwt/2 ;
        uppercurvfactor = uppercurvfactor/2 ; %%%%This implies (1/4) of the original value
        steps           = 3 ;
        stepsize        = 1/2 ;
        if length(knotvec_c)>=10
            pigrid_supp   = startvalgeneratorMultiStep(knotvec_c,pi_c,scalefactor,upperlinwt,uppercurvfactor,steps,stepsize,4,startvaldisplay) ;
        elseif length(knotvec_c)==9
            pigrid_supp   = startvalgeneratorMultiStep(knotvec_c,pi_c,scalefactor,upperlinwt,uppercurvfactor,steps,stepsize,3,startvaldisplay) ;
        elseif length(knotvec_c)==8
            pigrid_supp   = startvalgeneratorMultiStep(knotvec_c,pi_c,scalefactor,upperlinwt,uppercurvfactor,steps,stepsize,2,startvaldisplay) ;
        end
        %%%%Weed out any repeats:
        pigrid_supp = unique(pigrid_supp','rows','stable')' ;
        %%%%Weed out any vectors that violate feasibility:
        [eqindx,ineqindx,boundindx]  = feasiblecheck(pigrid_supp,Aeq,beq,Aineq,bineq,lb,ub,feastol) ;
        pigrid_supp(1,eqindx==0) = beq(1) ;
        pigrid_supp(2,eqindx==0) = beq(2) ;
        if sum(ineqindx==0|boundindx==0)>0
            pigrid_supp = pigrid_supp(:,ineqindx==1&boundindx==1) ;
        end
        ssrlog_supp   = nan(length(pigrid_supp(1,:)),1) ;
        disp(strcat('evaluating_',num2str(length(ssrlog_supp)),'_candidate parameter vectors.'))
        tic
        parfor ii=1:length(pigrid_supp(1,:))
            try
                pitemp = pigrid_supp(:,ii) ;
                ssrtemp = ObjFun_cMSM(pitemp,knotvec_c,student_id,Q,T,...  %%%
                    ThetaE,contract,Q80freq,pay,SimTimeCum,... %%%
                    Activeids,Marginalids1,Marginalids2,...
                    Tdom,EQgivenT_lb,EQgivenT_ub,cutoffs,... %%%
                    betaq1,betaq2,betaq3,knotvec_q,0,BSstate) ;
                ssrlog_supp(ii)  = ssrtemp ;
            catch
                ssrlog_supp(ii) = inf ;
            end
        end
        toc
        pigrid(:,length(pigrid(1,:))+1:length(pigrid(1,:))+length(pigrid_supp(1,:))) = pigrid_supp ;
        ssrlog(length(ssrlog)+1:length(ssrlog)+length(ssrlog_supp))                  = ssrlog_supp ;
        nextindx = find(ssrlog==min(ssrlog),1,'first') ;
        if length(nextindx)>1
            if max(abs(pigrid(:,nextindx(1))-pigrid(:,nextindx(2)))./pigrid(:,nextindx(1)))<10^-7
                pinext = pigrid(:,nextindx(1)) ;
            else
                warning('Two minimizers found where relative difference is less than 10^-7...')
                pinext = pigrid(:,nextindx(1)) ;
            end
        else
            pinext   = pigrid(:,nextindx) ;
        end
        format shortg
        CheckAltGuesses = [pi_c pinext 100*(pinext-pi_c)./pi_c; ssrlog(1) ssrlog(nextindx(1)) 100*(ssrlog(nextindx(1))-ssrlog(1))/ssrlog(1)]   %#ok
        pi_c = pinext ;
    end
    %%%%Now, we re-do the grid search once using the best guess from the previous round, a smaller
    %%%%stepsize, and smaller start values:
    disp('%%%%555555555555555555555555555555555555%%%%')
    disp('%%%% ADAPTIVE START POINT SEARCH, ROUND 5...')
    disp('%%%%555555555555555555555555555555555555%%%%')
    scalefactor     = scalefactor/2 ;
    upperlinwt      = upperlinwt/2 ;
    uppercurvfactor = uppercurvfactor/2 ; %%%%This implies (1/8) of the original value
    steps           = 2 ;
    stepsize        = 1/4 ;
    if length(knotvec_c)>=10
        pigrid_supp   = startvalgeneratorMultiStep(knotvec_c,pi_c,scalefactor,upperlinwt,uppercurvfactor,steps,stepsize,4,startvaldisplay) ;
    elseif length(knotvec_c)==9
        pigrid_supp   = startvalgeneratorMultiStep(knotvec_c,pi_c,scalefactor,upperlinwt,uppercurvfactor,steps,stepsize,3,startvaldisplay) ;
    elseif length(knotvec_c)==8
        pigrid_supp   = startvalgeneratorMultiStep(knotvec_c,pi_c,scalefactor,upperlinwt,uppercurvfactor,steps,stepsize,2,startvaldisplay) ;
    end
    %%%%Weed out any repeats:
    pigrid_supp = unique(pigrid_supp','rows','stable')' ;
    %%%%Weed out any vectors that violate feasibility:
    [eqindx,ineqindx,boundindx]  = feasiblecheck(pigrid_supp,Aeq,beq,Aineq,bineq,lb,ub,feastol) ;
    pigrid_supp(1,eqindx==0) = beq(1) ;
    pigrid_supp(2,eqindx==0) = beq(2) ;
    if sum(ineqindx==0|boundindx==0)>0
        pigrid_supp = pigrid_supp(:,ineqindx==1&boundindx==1) ;
    end
    ssrlog_supp   = nan(length(pigrid_supp(1,:)),1) ;
    disp(strcat('evaluating_',num2str(length(ssrlog_supp)),'_candidate parameter vectors.'))
    tic
    parfor ii=1:length(pigrid_supp(1,:))
        try
            pitemp = pigrid_supp(:,ii) ;
            ssrtemp = ObjFun_cMSM(pitemp,knotvec_c,student_id,Q,T,...  %%%
                                   ThetaE,contract,Q80freq,pay,SimTimeCum,... %%%
                                   Activeids,Marginalids1,Marginalids2,...
                                   Tdom,EQgivenT_lb,EQgivenT_ub,cutoffs,... %%%
                                   betaq1,betaq2,betaq3,knotvec_q,0,BSstate) ;
            ssrlog_supp(ii)  = ssrtemp ;
        catch
            ssrlog_supp(ii) = inf ;
        end
    end
    toc
    pigrid(:,length(pigrid(1,:))+1:length(pigrid(1,:))+length(pigrid_supp(1,:))) = pigrid_supp ;  
    ssrlog(length(ssrlog)+1:length(ssrlog)+length(ssrlog_supp))                  = ssrlog_supp ;
    [~,uniqueindx] = unique(pigrid','rows','stable') ;
    pigrid = pigrid(:,uniqueindx) ;
    ssrlog = ssrlog(uniqueindx) ;
    %%%%One last time: weed out any vectors that violate feasibility:
    [eqindx,ineqindx,boundindx]  = feasiblecheck(pigrid,Aeq,beq,Aineq,bineq,lb,ub,feastol) ;
    pigrid(1,eqindx==0) = beq(1) ;
    pigrid(2,eqindx==0) = beq(2) ;
    if sum(ineqindx==0|boundindx==0)>0
        pigrid = pigrid(:,ineqindx==1&boundindx==1) ;
        ssrlog = ssrlog(ineqindx==1&boundindx==1) ;
    end
    nextindx = find(ssrlog==min(ssrlog),1,'first') ;
    if length(nextindx)>1
        if max(abs(pigrid(:,nextindx(1))-pigrid(:,nextindx(2)))./pigrid(:,nextindx(1)))<10^-7
            pinext = pigrid(:,nextindx(1)) ;
        else
            warning('Two minimizers found where relative difference is less than 10^-7...')
            pinext = pigrid(:,nextindx(1)) ;
        end
    else
        pinext   = pigrid(:,nextindx) ;
    end
    format shortg
    CheckAltGuesses = [pi_c pinext 100*(pinext-pi_c)./pi_c; ssrlog(1) ssrlog(nextindx(1)) 100*(ssrlog(nextindx(1))-ssrlog(1))/ssrlog(1)]   %#ok
%     pi_c = pinext ;

    [~,keepindx] = mink(ssrlog,keepnum) ; 
    pigrid       = pigrid(:,keepindx) ;
    if nargout==2
        ssrlog   = ssrlog(keepindx) ;
    elseif nargout==3
        ssrlog   = ssrlog(keepindx) ;
        temp     = [pigrid0 pigrid] ;
        %%%%We started with nothing by pigrid0, and by appending pigrid0 to the final pigrid, the only way for one of the
        %%%%original columns to survive is if it appears twice in temp.  Therefore, length(uniqueindx)-keepnum has to equal
        %%%%the number of columns of pigrid0 that got replaced by a parameter vector with a better ssr value.  
        [~,uniqueindx] = unique(temp','rows') ; 
        survivorfrac   = (length(uniqueindx)-length(pigrid0(1,:)))/keepnum ;  
    end


end



